home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / LUDCMP.DEM < prev    next >
Text File  |  1991-04-29  |  3KB  |  101 lines
  1. PROGRAM d2r2 (input,output,dfile);
  2. (* driver for routine LUDCMP *)
  3. LABEL 10,99;
  4. CONST
  5.    np=20;
  6. TYPE
  7.    glnpbynp = ARRAY [1..np,1..np] OF real;
  8.    glnarray = ARRAY [1..np] OF real;
  9.    glindx = ARRAY [1..np] OF integer;
  10. VAR
  11.    j,k,l,m,n,dum : integer;
  12.    d : real;
  13.    a,xl,xu,x : glnpbynp;
  14.    indx,jndx : glindx;
  15.    dfile : text;
  16.  
  17. (*$I MODFILE.PAS *)
  18. (*$I LUDCMP.PAS *)
  19.  
  20. BEGIN
  21.    glopen(dfile,'matrx1.dat');
  22. 10:   readln(dfile);
  23.    readln(dfile);
  24.    readln(dfile,n,m);
  25.    readln(dfile);
  26.    FOR k := 1 to n DO BEGIN
  27.       FOR l := 1 to n-1 DO read(dfile,a[k,l]);
  28.       readln(dfile,a[k,n])
  29.    END;
  30.    readln(dfile);
  31.    FOR l := 1 to m DO BEGIN
  32.       FOR k := 1 to n-1 DO read(dfile,x[k,l]);
  33.       readln(dfile,x[k,l])
  34.    END;
  35. (* print out a-matrix for comparison with product of lower *)
  36. (* and upper decomposition matrices *)
  37.    writeln('original matrix:');
  38.    FOR k := 1 to n DO BEGIN
  39.       FOR l := 1 to n-1 DO write(a[k,l]:12:6);
  40.       writeln(a[k,n]:12:6)
  41.    END;
  42. (* perform the decomposition *)
  43.    ludcmp(a,n,np,indx,d);
  44. (* compose separately the lower and upper matrices *)
  45.    FOR k := 1 to n DO BEGIN
  46.       FOR l := 1 to n DO BEGIN
  47.          IF (l > k) THEN BEGIN
  48.             xu[k,l] := a[k,l];
  49.             xl[k,l] := 0.0
  50.          END ELSE IF (l < k) THEN BEGIN
  51.             xu[k,l] := 0.0;
  52.             xl[k,l] := a[k,l]
  53.          END ELSE BEGIN
  54.             xu[k,l] := a[k,l];
  55.             xl[k,l] := 1.0
  56.          END
  57.       END
  58.    END;
  59. (* compute product of lower and upper matrices for *)
  60. (* comparison with original matrix *)
  61.    FOR k := 1 to n DO BEGIN
  62.       jndx[k] := k;
  63.       FOR l := 1 to n DO BEGIN
  64.          x[k,l] := 0.0;
  65.          FOR j := 1 to n DO BEGIN
  66.             x[k,l] := x[k,l]+xl[k,j]*xu[j,l]
  67.          END
  68.       END
  69.    END;
  70.    writeln('product of lower and upper matrices (rows unscrambled):');
  71.    FOR k := 1 to n DO BEGIN
  72.       dum := jndx[indx[k]];
  73.       jndx[indx[k]] := jndx[k];
  74.       jndx[k] := dum
  75.    END;
  76.    FOR k := 1 to n DO BEGIN
  77.       FOR j := 1 to n DO BEGIN
  78.          IF (jndx[j] = k) THEN BEGIN
  79.             FOR l := 1 to n-1 DO write(x[j,l]:12:6);
  80.             writeln(x[j,n]:12:6)
  81.          END
  82.       END
  83.    END;
  84.    writeln('lower matrix of the decomposition:');
  85.    FOR k := 1 to n DO BEGIN
  86.       FOR l := 1 to n-1 DO write(xl[k,l]:12:6);
  87.       writeln(xl[k,n]:12:6)
  88.    END;
  89.    writeln('upper matrix of the decomposition:');
  90.    FOR k := 1 to n DO BEGIN
  91.       FOR l := 1 to n-1 DO write(xu[k,l]:12:6);
  92.       writeln(xu[k,n]:12:6)
  93.    END;
  94.    writeln('***********************************');
  95.    IF eof(dfile) THEN GOTO 99;
  96.    writeln('press return FOR next problem:');
  97.    readln;
  98.    GOTO 10;
  99. 99:   close(dfile)
  100. END.
  101.